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Abstract 

We derive and introduce anisotropic effective pair potentials to coarse-grain solutions of 
semiflexible rings polymers of various lengths. The system has been recently investigated by 
means of full monomer-resolved computer simulations, revealing a host of unusual features 
and structure formation, which, however, cannot be captured by a rotationally-averaged effec¬ 
tive pair potential between the rings’ centers of mass [M. Bernabei et al., Soft Matter 9, 1287 
(2013)]. Our new coarse-graining strategy is to picture each ring as a soft, penetrable disc. We 
demonstrate that for the short- and intermediate-length rings the new model is quite capable of 
capturing the Physics in a quantitative fashion, whereas for the largest rings, which resemble 
flexible ones, it fails at high densities. Our work opens the way for the physical justification of 
general, anisotropic penetrable interaction potentials. 
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1 Introduction 


By the simple process of joining the ends of a linear polymer chain, one obtains a ring polymer 
(RP) 1 . While the architecture of ring polymers is very simple, they differ in many interesting 
ways from their linear counterparts and are the subject of active research in Physics, Biology, 
Chemistry and even pure Mathematics. One interesting consequence is that for a dynamics that 
disallows strand crossing, there are different classes of configurations of a RP, which can never 
transform into each other. These are referred to as the topology classes or knot types of a RP. Knot 
theory is a fascinating and active branch of mathematics with many open problems concerning 
the enumeration and classification of knots®. A RP that has the topology of a circle is called an 
unknotted RP. 

Unlike RPs, linear polymer chains can strictly speaking never be knotted, as every configura¬ 
tion of a linear polymer chain can be continuously transformed to a straight line, without the need 
for strand crossings. While this is a fundamental difference between linear polymer chains and 
RPs, it is nevertheless possible to extend the concept of knots to physical knots on linear poly¬ 
mer chains 3 . There are many works dealing with the properties of these physical knots on linear 
polymer chains^ 1 0, in particular due to their relevance in Biophysics, where they are for instance 
found in DNA 12 13 and can have significant effects on key processes 1416 . 

The topological constraint of a ring polymer has important consequences for its physical be¬ 
havior. It took the work of many authors 17 24 to establish that the diameter of gyration D g of an 
isolated ideal RP with fixed topology scales as (D 2 ) ~ N 2v , where v ~ 0.588 is the Flory exponent, 
which also describes the scaling behavior of the radius of gyration of self-avoiding linear chains 25 . 
This is true for all knot types, even for an ideal unknotted RP (i.e., without monomer excluded vol¬ 
ume interactions, just keeping the topological constraints). An ideal linear polymer chain, on the 
other hand, remains more compact and exhibits the scaling law of a non-self-avoiding random walk 
(D 2 ) ~ N. Another important difference lies in the effective potential between the centers of mass 
of RPs. While the effective potential vanishes between infinitely thin linear polymer chains and 
also between an infinitely thin linear polymer chain and a RP, there remains a nonzero repulsive 
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contribution between cyclic polymers with fixed topology 26 27 . Here one usually speaks of a topo¬ 
logical potential. Furthermore, it was shown that the effective potential between two moderately 
sized RPs increases with the knot complexity of the rings 27 . 

Also for concentrated systems, the topology of polymer chains plays an important role. The 
scaling of linear polymer chains in the melt is the one of a non-self-avoiding random walk (D;) ~ 
A. Simulations of dense systems of RPs 28 29 on the other hand showed that while short chains 
also exhibit a Gaussian scaling behavior (D 2 ) ~ A, long chains are compact and thus scale as 
(D 2 ) ~ A 2 / 3 . In between there is a broad crossover region, where a (D 2 ) ~ A 4 / 5 scaling provides 
a good description of the data. For the dynamics, it is expected that concatenations of ring poly¬ 
mers can have a significant effect, as they are permanent, in contrast to the entanglement of linear 
polymer chains. However, this implies that those concatenations are there in the first place, i.e., 
from the very synthesis of the sample on. Even in the absence of concatenations, there are impor¬ 
tant differences in the dynamics of RPs in the melt with respect to their linear counterparts. For 
instance recent experiments 30 and simulations 31 revealed a power-law stress relaxation instead of 
the rubbery plateau found for linear chains. 

For the large intermediate density domain between dilute solutions and melts, there are rela¬ 
tively few theoretical results despite the practical relevance of this regime for instance in the field 
of biophysics, where the topological interactions between chromatin loops plays a crucial role in 
the creation of chromosome territories 26 32 34 . A fruitful and modem approach for the economic 
description and simulation of macromolecules in this regime is the method of coarse-graining. 
The idea behind this method is to bridge the time and length scales in the system by describing the 
macromolecules via an effective model with a reduced set of suitably-chosen effective degrees of 
freedom (d.o.f.). The microscopic information of the monomer-resolved model is underlying the 
effective model, as it determines the form of the effective potential, which describes the interaction 
between the macromolecules. The advantage of this method is not only that every timestep in a 
simulation requires less computational effort due to the simplified representation, but also that one 
can often choose a much larger timestep in a simulation of the coarse grained model, as the d.o.f. 
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that remain in the coarse-grained model change much slower in time than their counterparts in the 
monomer-resolved model 35 36 . 

The method of coarse-graining is well-established and has for instance found successful appli¬ 
cation for polymer chains 37 39 i star polymers 40 ' 43 ^, star-shaped poly electrolytes 44 45 , dendrimers 4648 , 
and block copolymers 49 51 . The identification of the relevant degrees of freedom is an essential 
part in the design of an effective model. One often uses isotropic effective models, where the 
macromolecules consisting of many individual monomers are reduced to their center of mass. For 
the semiflexible ring polymers such a model has already been investigated in ref. 52 . Clustering 
was observed in monomer-resolved simulations of semiflexible ring polymers, as well as in the 
corresponding isotropic effective model. However, it was also shown that the monomer-resolved 
system shows anisotropic features that can not be accounted for in the isotropic effective model. 
Also the correlation functions stemming from the isotropic effective model are markedly different 
from the microscopically derived ones. Anisotropy is particularly strong for rings with high bend¬ 
ing stiffness or few monomers, as they have a strong tendency to orient with respect to other rings 
in their proximity. This motivates us to introduce an anisotropic effective model for the description 
of semiflexible ring polymers in this article. In this model, we will define the effective particles 
as soft disc-like molecules which are described not only by their center of mass but also by the 
direction in which their faces are oriented. An anisotropic effective model was already used suc¬ 
cessfully for the description of hard disc-like macromolecules 53 , but to the best of our knowledge 
this approach was up to now never applied to penetrable macromolecules, where the centers of 
mass of the macroparticles can coincide. Penetrable particles are particularly interesting as they 
allow for clustering, which often leads to a rich phase behavior. For instance, point-particles in¬ 
teracting with a certain class of ultrasoft potentials form so-called cluster crystals® 4 56 . Unlike 
in an ordinary crystal multiple soft particles can sit on top of each other at the same lattice site 
in a cluster crystal. Another peculiar feature of this state of matter is that by compressing it one 
only changes the occupation number of particles per lattice site, while the lattice constant remains 
invariant. Monomer-resolved simulations of semiflexible ring polymers on the other hand show 
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the formation of the cluster glass phase 57 58 , which is an arrested-state that also contains some of 
the features found in the cluster crystal phase. In both cases the overall structure of the system is 
frozen, while individual particles can hop between the lattice sites of the cluster crystal or the stacks 
found in the cluster glass. Elongated dendrimers, which unlike hard rod-like particles exhibit local 
antinematic order 5 ^, are another interesting example for a a system of penetrable particles, which 
behaves distinctively different to its solid counterpart. By creating an anisotropic effective model 
for the semiflexible ring polymers, we aim at a model that is still computationally cheap, and that 
improves the description obtained by the isotropic model, especially for the case of high densi¬ 
ties. In addition, the analysis of the interactions in the anisotropic effective potential allows us to 
get a better understanding for the interaction between the anisotropic, penetrable nanoparticles in 
systems of semiflexible RPs. 

The remainder of this article is structured as follows. In Section [2] we first present the Hamil¬ 
tonian of the monomer-resolved model which we use for the description of semi-flexible RPs and 
then introduce an anisotropic effective model for such a system. In Section [A] we give more details 
about the derivation of the effective interactions for such a model. We carried out Molecular Dy¬ 
namics simulations of the monomer-resolved model and Monte Carlo simulations of the effective 
models; details concerning these simulations are given in Section [3] In Section [4] we present the 
anisotropic effective potential and discuss its features. Results of Monte Carlo simulations with 
this potential, which show that the inclusion of anisotropy in the effective model can significantly 
improve the agreement with the monomer-resolved model, are presented in Section [5j whereas in 
Section[6]we briefly discuss the effects of truncation of the expansion of the potential on the quality 
of the results. Conclusions are given in Section [7} In the Appendix, we explain the expansion of 
the anisotropic pair-correlation function of a system of two RPs, which contains all the information 
for calculating the effective potential, as a sum of suitably chosen basis functions. 
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2 Anisotropic Effective Model 


2.1 Monomer-resolved Model for Semiflexible Ring Polymers 

The derivation of the anisotropic effective model is based on a microscopic model of semiflexible 
ring polymers, each consisting of N monomers. They are described with the bead-spring model 
by Kremer and Grest 60 and an additional rigidity term. Thus, any two-monomers interact via the 
truncated and shifted Lennard-Jones potential 


V LJ (r) 




if r < 2 1 / 6 cj; 


0 


if r > 2 1 / 6 o\ 


( 1 ) 


This potential is purely repulsive, accounting then for monomer excluded volume interactions. 
Bonded monomers also interact through a finitely extensible non-linear elastic potential (FENE) 


VfeneM = —^ In 



( 2 ) 


Rigidity is introduced via the bending potential 


V b end(e) = K(l-COS0) 2 , (3) 

where 0 is the angle between two consecutive bond vectors Q The potential Vbend vanishes for 
0=0, when the polymer chain does not bend at the respective angle. We choose £ = k B T, k = 
30 k B T / a 2 , Rq = 1.5a and K = 30 k B T, where k B is the Boltzmann constant and T the temperature. 
These are precisely the parameters employed in the simulation study of ref. 52 . The corresponding 
dynamics does not allow for chain crossings and thus topology is preserved. 

In ref. 58 the characteristic ratio 61 of this polymer model was estimated by carrying out sim¬ 
ulations of isolated linear chains. Excluded volume interactions were switched off except for 

'Note that this is not the Kratky-Porod model (linear in the cosine). We expect the same qualitative results for 
Kratky-Porod rings with the same N’s and persistence lengths as the model simulated here" 
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mutually connected monomers, in order to obtain long-range Gaussian statistics. CA was obtained 
by analyzing the long-5 limit of the ratio (R 2 (s))/s(b 2 ), where R(s) is the distance between two 
monomers i,j with s = |z — j\, and b is the bond length ((b 2 ) = 0.94). The authors reported a 
value of Coo ~ 15, which is typical for stiff polymers® 1 ! We can give an estimate for the per¬ 
sistence length of the model by mapping it to the freely rotating chain model using the relation 
cos0 = ^°°~| ~ 0.875, where 6 is the bending angle of the freely rotating chain model 61 . The 
persistence length is then obtained as s p b — — &/ln(cos0) ~ 7.3. We carried out simulations of 
ring polymers with N — 20, 50, and 100 monomers, which have the contour to persistence length 
ratio N/sp ~ 2.7, 6.7 and 13.3 respectively. 


2.2 The Anisotropic Effective Model 

In earlier work 52 , Bemabei et al. carried out extensive monomer-resolved simulations of stiff ring 
polymers to obtain the structure of concentrated solutions of the same. In an attempt to coarse- 
grain the system in the simplest possible way, they also derived and employed an isotropic effective 
potential for their effective description, reducing thereby stiff RP’s into point-like effective parti¬ 
cles, namely their centers of mass. At this level of approximation, the effective particles possess 
no other, internal (spin-like) degrees of freedom and thus the effective interaction is isotropic. The 
effective potential between these macroparticles was defined by calculating the pair correlation 
function g lso (r) in an infinitely dilute system and using 


j3Vg(r) = ~In g iso (r) 


(4) 


to define the effective interaction potential between these point particles, where /3 = 1/ (k^T). In 
the infinitely dilute case, the distribution of the centers of mass of the ring polymers in equilibrium 
is identical to the distribution of the point particles in the effective model. At higher densities, 
however, it turns out that multi-particle terms in the effective potential are necessary to obtain 
the correct equilibrium distribution of the centers of mass in the effective model. One of the 
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reasons due to which multi-particle interactions become important is that two ring polymers that 
are sufficiently close and stiff, will prefer to align parallel to each other. A third ring polymer 
interacting with those two will not see them as two independent rings but as a system of two 
rings that are correlated. When using the potential V e ff(r) calculated in Eq. ([4]) one assumes that 
the free energy penalty of one ring with respect to a second is independent of the presence of 
another polymer in the vicinity of the second. If the ring polymers preferentially align parallel, 
this assumption is clearly violated and one has to correct the effective potential by introducing 
multi-particle terms. Bemabei et al. showed that already at moderate densities one can encounter 
strong correlations of the orientations of the semiflexible ring polymers, in particular if the chains 
contain only few monomers (e.g., N = 20) 52 . Therefore, a more accurate coarse-graining which 
takes into account the ring anisotropy is called for. 

The easiest way to incorporate the correlations of the orientations of rings is to introduce them 
via additional degrees of freedom (d.o.f.) in the effective description. This is precisely what we 
do in this article. For this purpose, we need to first come up with a suitable definition for the 
orientation of a RP. To this end, we make use of the gyration tensor 


*^ 0!/3 


1 

N 


N 

E 

i=l 


— > 77/ r, 


(0 JO 


(5) 


where rjy (a = x,y,z, Cartesian components) denotes the position of the i-th monomer with respect 
to the center of mass of the ring to which this monomer belongs. The eigenvectors of S a p are 
the principal axes of an ellipsoid that approximates the shape of the macromolecule: If a RP is 
flat, which means that all its monomers lie in one plane, the ellipsoid has one zero eigenvalue 
with a corresponding eigenvector that is perpendicular to that plane. Also in the more general 
case, where the monomers do not all lie in the same plane, we define the normalized eigenvector 
corresponding to the smallest eigenvalue of S aj g as the direction vector d of the RP. Note that d and 
—d are equivalent for reasons of symmetry. The ring polymers in the anisotropic effective model 
we propose are described via the position vectors of their centers of mass, RW, and their direction 
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vector d ( T Henceforth, we describe the stiff rings as soft circular discs and ignore differences in 
the other 2 eigenvectors of the gyration tensor S a p . This choice is motivated by the limit of infinite 
bending stiffness, where the rings assume flat and precisely circular conformations. Our model, 
therefore, amounts to the minimal anisotropic extension of the spherically symmetric effective 
interaction between the centers of mass. We emphasize, however, that there is no a priori guarantee 
that this will be an improvement over the isotropic model at finite densities and in particular at high 
concentrations: this depends on the degree in which the RP’s at high concentration maintain their 
anisotropic shape and properties encoded in the high-dilution limit in which the anisotropic pair 
potential is derived. Accordingly, the introduction of such a potential is not a straightforward part 
of a systematic strategy of introducing more and more detail into the effective description of the 
system. 

In order to determine the anisotropic effective potential V e ff, we carried out monomer-resolved 
simulations of two ring polymers inside a large simulation box. The effective pair potential is then 
defined such that it exactly reproduces the correlation functions of the effective degrees of freedom 
in this infinitely dilute, monomer-resolved simulation. In the effective model, two ring polymers 
are described by a total of 10 degrees of freedom, three for each center of mass and two for each 
direction vector of each ring polymer. However, due to translation, rotation and mirror symmetry 
the distinct configurations (those that cannot be related by symmetry transformations) of a system 
with two effective particles are reduced and can be specified by 4 parameters only. 

A convenient choice for these variables is illustrated in figure |T] and reads as follows: 

r = | r |; 
cos 0i = d ( 1 ^-r; 
cos 0 t = d ( 2 ^-r; 

9 ' arccos (|d)VAJ’ ( ) 

where r = R^ 2 )—R (1 ) is the connection vector between the centers of mass of the two rings, r = r/r 


9 



the unit vector in the direction of r and the component of the director d ''- 1 perpendicular to r; 
0 < (p < K denotes the angle between vectors and d “ . By selecting the appropriate sign of 
d !,) we can always choose cos 0, to lie in the interval [0,1]. 



Figure 1: Illustration of the effective variables r = |r|, 0, and (p with which relative configurations 
of ring polymers are described. 

Let us define the ideal case as the system where the effective particles do not interact and thus 
every orientation and position of the effective particles occurs with equal probability, independently 
of the configuration of the other effective particle. With the effective coordinates defined in (|6]) the 
probability density in the ideal system, (r, cos 0], cos 02 , <p) is proportional to r 2 and constant in 
both cos 0, as well as in (p. This simple behavior of Pid(r, cos 0i, cos 02, <f>), makes ([6]) a particularly 
convenient choice of the effective coordinates. In the simulations with two ring polymers, we 
obtain a probability density P (r,cos 0| ,cos Qi-(p) that is different from the ideal distribution Pjj. 
We define as a generalised version of the radial distribution function the anisotropic pair correlation 
function as 


g (r, cos 0i, cos 0 2 , <p) 


P(r,cos 0i, cos 02, <p) 
Pid (r, COS 01 , COS 02 , (p) 


(7) 


Thus, the quantity g (r,cos 0i,cos 02, (p) describes the factor by which configurations in the effec¬ 
tive anisotropic model have to be enhanced or suppressed with respect to the ideal case, in order 
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to obtain a distribution for the effective d.o.f. that is identical to the distribution obtained in a 
monomer-resolved simulation in the infinitely dilute case. As in the isotropic case 0 the relation 
to the associated effective potential reads 

j6V e ff (r,cos 0i, cos 02 , <p) = - In [g(r, cos 0i, cos 0 2 , <?>)]■ (8) 

From the anisotropic effective potential, we can deduce the isotropic pair-correlation function via 

g lso (r) — “ / dcos0i / dcos0 2 / d^ g(r,cos0i,cos02,<p). (9) 

Ti Jo Jo Jo 

The associated isotropic effective potential is then given by 0. 

The effective potential between two identical ring polymers remains invariant if we swap the 
orientations of their respective director with respect to the connection vector, i.e., if we swap the 
values of the polar angles 0i and 02 of the two rings: 

V e ff (r,COS0i,COS02,<p) = Veff(r,COS02,COS0i,<p). (10) 

This symmetry is violated for the effective potential between bidisperse ring polymers, e.g. for 
ring polymers with a different number of monomers. However, apart from this symmetry, there 
would be no differences in the procedure of calculating the effective potential between different 
types of ring polymers. 

3 Simulation Details 

3.1 Derivation of the effective interaction 

To determine the effective potential, we carry out constant NVT molecular dynamics simulations 
with two ring polymers. We simulate rings of N — 20, 50 and 100 monomers. For these simulations 
we use the LAMMPS simulation package 62 . The polymer rings are placed in a simulation box, 
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which is large enough to prevent multiple interactions via the periodic boundary conditions. The 
temperature in the simulation is maintained by the use of Langevin-Dynamics. The corresponding 
equations of motion read as®: 


mfi(t) = F i{t) - ymviit) + rift). (11) 

Here, r, is the position of the i-th monomer, m its mass and F,(/) the deterministic force acting on 
it, which includes the microscopic forces originating from potentials (|T[j3]) and the force originating 
from a biasing potential. The bias potential Vbi as = { r ~ r j) kj/2 introduces a harmonic spring 
with spring constant kj between the centers of mass of the ring polymers. The spring is relaxed for 
r — rj. We carried out simulations for different values of rj, starting at r f = 0 and increasing it up 
to some maximum value rc in steps of a/ 2. VoxN — 20, 50, 100, rc was chosen as 10(7, 20c, 30c 
respectively. These values for r c are much bigger than the infinite-dilution diameters of gyration 
(D g o = 5.9(7, 13(7, 21.5(7 for N = 20, 50, 100 respectively). For kj we chose the values 2.5e/a 2 
and e/o 2 for all ring sizes and for N = 20 we also carried out simulations with kj — 5e/o 2 . The 
quantity 7] ft) is a random force, with (rift)) = 0, which is related to the friction coefficient y 
by the fluctuation dissipation relation it')) = 2 yinksT8ij8 a p8(t — t'), a and /3 denoting 

Cartesian components. Our unit of time is set by to — (ma 2 /e) 1 / 2 and the friction coefficient y 
is chosen as 1 /to. We integrate the equations of motion with a timestep of At = 10 3 /o, and use 
2 x 10 8 timesteps for equilibrating the system and collect data during another 2 x 10 9 timesteps. 

We sample histograms ( Q ), where Q refers to a bin in the 4D space of the effective coordi¬ 

nates. P U \Q) gives the probability for a state in the j-th simulation to have effective coordinates 
in bin Q. The histograms have 128 bins in r and 16 bins in cos 66, cos 02 and (p direction. As 
discussed in appendix |A| we use the Self-Consistent Histogram Method by Ferrenberg and Swend- 
sen 6465 to combine the P ij HQ) histograms for simulations with identical ring sizes but different 
biasing potentials to arrive at an estimate for P(Q) in the unbiased system. 
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3.2 Many-body effective fluid 


Using the anisotropic effective potential we carry out standard Metropolis Monte Carlo (MC) sim¬ 
ulations for the anisotropic effective model. The values of the anisotropic effective potential have 
been calculated on a discrete grid in the (r,cos Q \,cos 02, (p) space and we use linear interpolation 
to estimate the values of exp(— /3 V e ff) in between the grid points. Having in mind a comparison 
with both the monomer-resolved simulation results and the isotropic effective potential of ref. 5 ^, 
we choose the same number of particles and effective densities that were used in those simulations. 
For rings with N = 20, 50 and 100 monomers we simulate systems of n = 2400, 1600 and 1200 
rings, respectively, varying in each case accordingly the cubic box size L as to achieve the desired 
density p =n/L 3 . As the effective potential V e ff is bounded, a random distribution of the particles 
in the simulation box can be used as initial condition. We have implemented two types of MC 
moves: the first one translates a randomly chosen particle in a random direction, and the second 
randomly rotates the particle’s director by some angle. The distance by which the particles are 
displaced and the angle by which they are rotated is randomly selected in an interval starting at 
0 and going up to some maximum value. For both moves, this maximum value of the interval is 
chosen such that the acceptance ratio is approximately 15%. We use 9 x 10 6 MC moves to equi¬ 
librate the system. During this equilibration period the individual soft particles diffuse to several 
times their own diameter. Afterwards, during 15 x 10 6 MC moves equilibrium configurations are 
generated. We store the configurations every 20 x 10 3 moves and use them to compute the physical 
observables that are presented in the section [5} 

The gain in computational efficiency for the simulations in going from a monomer resolved to 
a coarse grained simulation is considerably. The relevant quantity to consider here, would be the 
velocity through phase-space. This can conveniently be characterized by means of the mean square 
displacement of the rings per unit of CPU time. If one ignores the detailed implementation aspects, 
the CPU time spend on a single sweep over all monomers in the former and a run over all effective 
ring particles in the latter, which strictly speaking depend on both the number of monomers N 
per ring and the overall density of rings, are of similar magnitude. However, the diffusion per 
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sweep in the coarse grained simulation is significantly larger than that for the monomer resolved 
simulation, i.e., for the case of TV = 50 and p* —20 this results in a factor of approximately 10 4 . The 
reason for this dramatic improvement is two-fold. First of all the translation/rotation of an effective 
ring corresponds to a much more time-consuming collective movement of the constituents. The 
second even more important contribution arises from the steric interaction that are present in the 
monomer resolved simulations and prevent the unphysical crossing of chain segments. In the 
coarse grained simulations such a restriction is absent, i.e., the effective rings are penetrable and 
can move apparently through each other. On this level of description this effect is not an unphysical 
process, but should be interpreted as a short-cut connecting initial and final configurations that 
are connected by a much more time-consuming and physically realizable pathway of folding and 
collective monomer movements. 


4 The Anisotropic Effective Potential 

We commence by recalculating the isotropic , i.e., angularly-average effective pair potential V^“(r) 
between the stiff rings, as a way of comparison with the previously derived results in refill Results 
are summarized in figure [ 2 ] reproducing indeed the previously derived ones 52 . The potential for 
r — 0 is finite, as the rings are allowed to overlap. It also features a local minimum there, whereas 
its maximum is located at r ^ 0.25D g o for all ring types investigated. Here, D g o is the average 
diameter of gyration of a free ring polymer. The height of the potential barrier at small distances 
of r decreases by increasing the number of monomers N on the ring polymers. 

Going over to the anisotropic effective model, we proceed with computing the aforementioned 
pair correlation function g(r, cos 0], cos 02, (p ) for a system of two ring polymers. As the latter de¬ 
pends on 4 effective coordinates and is therefore difficult to visualize, we first introduce a reduced 
pair correlation function g(r. d 1:11 • d' 2 - ) ), which expresses the relative joint probability density of 
observing the two ring polymers at a distance r and with the directors mutually oriented at the 
value given by their scalar product, d* 1; ■ d :2) , over the same quantity for noninteractive rings. 
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Figure 2: The effective, center-of-mass pair potential in the isotropic effective model for rings with 
different numbers of monomers N. The center-of-mass separation r is scaled with D a o, the average 
diameter of gyration of a free ring polymer. The solid line shows the angularly-averaged effective 
pair potential, while the dashed lines are results of refj 52 . 
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r /° E 0 

(c) N = 100 


Figure 3: Infinite-dilution limit of the quantity G(r, d 1 j ■ d i2j ), which quantifies the distribution of 
the scalar product between directors for different values of r. We visualize this distribution for the 
ring sizes N — 20, N = 50 and N — 100. 
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Moreover, we introduce a reduced version of this function, G(r,d (1) ■ d 2J ), by dividing over its 
isotropic counterpart, i.e., 


G(r,d (1) -d (2) ) 


g(r,d (1) -d (2 ^) 

g 180 M 


( 12 ) 


Results are summarized in figure [3] One can see that the angular distribution of the directors 
changes significantly for r « 0.25D g o, which is approximately the position of the maximum of 
V-(r). For r < 0.25D a o the angle between the directors is biased towards %/ 2, while for 0.25D g o < 
r < D g o they prefer to align parallel with respect to each other. The position of the maximum of 
Vg S f °(r) coincides approximately with the distance r where interpenetrated configurations of the 
rings become sub-dominant and where they are more likely to align parallel to each other. The 
transition between these two domains is particularly steep for N — 20 and becomes smoother for 
rings with a larger number of monomers. When the rings interpenetrate each other, the distribu¬ 
tion of angles between the directors is rather wide, while it gets narrow after the transition where 
the bias towards parallel alignments of the rings is very strong in particular for the smallest rings 
with N — 20. It is readily visible from figure |3]that anisotropy is particularly important for smaller 
rings. For r > D g o, the distribution of the angle between the directors becomes flat, as the rings 
are then well separated and hence do not interact. Note that by definition, Eq. ( [T2| ), the quantity 
G(r,d^ • d f2) ) is a normalized probability distribution for fixed r, and in figure [ 3 ] it is therefore 
meaningless to compare the plotted function at different r values. 

The relative orientation between the vectors r, d^ and d' :2i is of course not completely deter¬ 
mined by the scalar product d 11 ^ • d' 2 -*; the function G(r. d' : 1 ' • d !2 *) contains less information than 
the full correlation function g(r, cos B \, cos 02 , <f>). In particular, when the directors are parallel, i.e., 
d ( 1 ). d < 2 ) = 1, the angle between the connection vector r and the directors d (1 ) and d (2 ) is still arbi¬ 
trary. We denote a configuration with d (1 ) || d* 2 * || r as-and a configuration with d^ 1 * || d (2 ^ _L r 

as ||. From the reduced pair correlation function G(r,d*^ -d (2 )) alone, we cannot say which of 
these two configurations is more probable, as the scalar product d ( 1 ' • d ! 2 * is identical to 1 in both 
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Figure 4: The effective potential for three different, fixed configurations of the directors and the 
connecting vector. As a comparison we also plot the pair potential in the isotropic effective model, 
(a) N = 20; (b) N = 50; (c) N = 100. The effective potentials are shown only for r values for which 
we have relatively good statistics. We also show a sketch of the 11,-and |— configurations. 
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cases. Using the full anisotropic pair correlation function g(r, cos 6 \, cos 0 2 , 9 ) we can compare the 
corresponding effective potentials: 

PV —(r) = — In [g(r,cos Q\ = l,cos0 2 = 1 , 9 )] 

P v \\( r ) = -ln[g(r,cos0i =0, cos 0 2 = 0,9 = 0 )]. (13) 

For the-case the value of the 9 coordinate is immaterial. However, due to the finite bin 

size of the grid on which we have calculated g(r,cos 0 i,cos 0 2 , 9 ) the choice of 9 makes a small 

difference, even for V _(r). We compute V _(r) from the average of g(r, cos 0i = 1, cos 0 2 = 1, 9 ) 

in 9 . 

In figure |4j we see that Vj| (r) increases significantly when r approaches D g o, while V _stays 

close to 0 until much smaller distances r. We can understand this results if we imagine the rings as 
discs with diameter D g o- In the 11 configuration, the rings lie in the same plane and will therefore 
start to overlap as soon as r < D g o- Since the rings are not perfect circles and their shape fluctuates, 
they can feel each other also for distances r which are slightly larger than D g o- In the-config¬ 

uration two discs overlap only if the distance between their centers of mass is smaller than their 
thickness. These results tell us that the peak in the reduced pair correlation function g(r, d 11: ■ d' :2 ' ) ) 

for r ~ 0.25D a o and d (1) • d (2 ^ « 1 is mostly due to-like configurations. However, as soon 

as the rings can overlap in-type configurations the effective potential increases very fast for 

smaller r and we come to a regime where other configurations of the directors are more favorable. 
As a comparison we also consider a configuration with d (1 ^ _L d (2 ^ || r, which we denote by | —. 
The corresponding effective potential is given by 

1 8 Vj_(r) = -ln[g(r,cos0i =O,cos0 2 = 1 , 9 )]. (14) 

As in the-case, the value of the 9 coordinate is irrelevant for calculating Vj_(r), which we 

compute from the arithmetic mean of g(r,cos 0\ = 0 ,cos 0 2 = 1 , 9 ) in 9 . 

For small distances r one ring interpenetrates the other in microscopic configurations of type 
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|—. While Vj_(r) starts to increase at larger r values than V _(r), the increase is slower and con¬ 

verges to a constant for r —* 0. This is intuitive to understand since it requires only a finite amount 
of bending energy to deform two rings such that one can fit into the other. The required bending 

energy is smaller if the rings are larger. The | — becomes dominant over the-configuration at 

an r value below the threshold r ~ 0.25D g o. This is also the r value at which we find the transi¬ 
tion in the reduced pair correlation function g(r,d ( ^ • d i2j ) between a regime where configurations 

with parallel directors, as in-, are preferred, to a regime where they are suppressed and other 

configurations like |— become dominant. 

In the Appendix we explain how one can expand the angular part of g(r,cos Q\ .cos O 2 , (p) into 
a series of suitably chosen basis functions f/j,/ 9 , w (cos0i,cos 02 , <P)- The expression for this ex¬ 


pansion is given in Eq. (27) and the corresponding coefficients c i u i i m {r) can be determined by 


calculating particular ensemble averages as shown in Eq. (30). We plot these coefficients c i u i 2 tn {r) 
for /i ,/2 < 2 in figure [5] From the fast change of c/ l 5 / 2 , m (r) for r ~ 0.25D g o one can once more 
see the transition between two regimes for r in which the distribution of the directors of the ring 
polymers is very different. We can again see that this transition is smoother for larger rings. The 
magnitude of coefficients c/ lj /, jW (r)/co,o,o(^) with (h,h) 7 ^ ( 0 , 0 ) tells us about the significance of 
the corresponding anisotropy in g(r,cos 0i,cos O 2 , <p). Anisotropy is more important for smaller 
rings and becomes more pronounced after the transition at r ^ 0.25D a o, where the rings prefer 
parallel configurations. 


5 Monte Carlo Simulations of the Anisotropic Effective Model 

We carried out Monte Carlo simulations of systems of effective particles described by the anisotropic 
effective model for different ring sizes N and various densities. We define the reduced density in 
our simulation as p* = «D 3 0 /L 3 , where n is the number of rings in the sample. In order to assess 
the quality of the anisotropic effective model, we compare our results to results of full monomer- 
resolved simulations from refJ 52 and the results of simulations using the isotropic effective model. 
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Figure 5: The first coefficients in the expansion of g(r,cos 6 \.cos 62 , (p ) divided by the coefficent 
c 0 ,o,o(>). (a) N = 20; (b) N = 50; (c) N = 100. 
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Figure 6: The pair correlation function g(r) at low reduced densities p* for a simulation of many 
ring polymers in the full monomer-resolved simulation (symbols), the anisotropic effective model 
(solid line) and the isotropic effective model (dashed line). 
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As we can see in figure [6] for all choices of the number of monomers N the effective models are 
in good agreement with the full monomer-resolved simulations at low densities p*. This is an 
important consistency check for the effective models, in which the interactions have been chosen 
such that the distribution of the effective degrees of freedom agrees with their distribution in the 
full monomer-resolved simulations, in the limit of small densities. 



Figure 7: The pair correlation function g(r), at high densities, for a simulation of many ring poly¬ 
mers with N — 20 monomers in the full monomer-resolved simulation (symbols), the anisotropic 
effective model (solid line) and the isotropic effective model (dashed line). 

In figure [7} we present results for the smallest rings with N — 20 monomers at higher densities. 
There is a dramatic improvement of the accuracy as one compares the isotropic with the anisotropic 
model. While the former fails for p* >2 the anisotropic effective model works up to p* = 5 and 
even gives a semi-quantitatively correct description of the system at p* = 5.97. At the highest 
densities, we see the development of a peak in g(r) at r = 0.3D g o. This peak in the pair correlation 
function is associated with the emergence of stacks of parallel rings and its position describes 
the typical distance of rings in these stacks^ 2 ! Interestingly, the isotropic effective potential has 
a maximum for r ~ 0.25D g o which is close to the typical distance of the rings in the stacks and 
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one could wonder why the rings prefer to align at a distance which seems to have a very high free 
energy penalty according to the isotropic effective potential. The answer to this apparent paradox 
lies in the strongly peaked nature of the anisotropic effective interaction, which we could observe in 
figures [3] [4] and [5] While the average configuration of the angular degrees of freedom at distances 
r ~ 0.3D g o has a high free energy penalty, a certain class of configurations, where the directors of 
the rings are almost parallel, is much more favorable. Obviously, stacking can not be observed in 
the isotropic effective model, where particles possess no directional degrees of freedom. 



Figure 8: (n r ) in a simulation of many ring polymers with N — 20 monomers in the full monomer- 
resolved simulation (symbols), the anisotropic effective model (solid line) and the isotropic effec¬ 
tive model (dashed line). 


As a further characteristic of the short-range coordination of the rings, we consider the average 
number (n r ) of neighbors within a distance r from the center of mass of a randomly chosen ring. 
This is expressed as 


(n r )=4np / dxx 2 g(x). 

Jo 


( 15 ) 
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For the rings with N — 20 monomers we present results for (n r ) in figure [8] Once more, the 
good agreement between the full monomer-resolved and the anisotropic effective model, even at 
the highest densities investigated, is confirmed: small differences appear only for 0.25D g o < r < 
0.4D„o- Evidently, (n r ) does not contain more information than the g(r) plot in figure [7j but it 
nevertheless clarifies the meaning of the disagreement between the g(r) curves in the monomer- 
resolved and the anisotropic effective model. At the highest densities in the full simulation, the 
centers of mass move a bit closer to each other than they do in the anisotropic effective simulation. 
This manifests itself as a shift of the peaks in the g(r) curves. The difference in the height of 
the peaks is partly a consequence of the shift, since a peak in g(r) has to be higher at smaller 
distances if it amounts to the same amount of average neighbours as a peak at a larger distance r. 
The fact that the (n r ) curves for the anisotropic effective and the full simulation in figure [8] agree 
for r > 0.4D a o shows us that the peaks in the g(r) curves indeed correspond to the same amount 
of average nearby particles that are simply accumulated at slightly different distances. 

We proceed now with the longer rings, N — 50. As can be seen by the pair correlation curves in 
figure [9j a), also in this case the inclusion of anisotropy improves the agreement with the monomer- 
resolved simulations significantly for densities p* from 2.3 to 10.2. In figure |9}b) we present g(r) 
for higher p*. In the full monomer-resolved simulations we can see a peak emerging in the pair 
correlation function g(r) at r = 0 on increasing the density. As described in ref. 52 58 the monomer- 
resolved system forms stacks of quasi-parallel oblate rings that are fully penetrated by bundles of 
elongated rings. In this phase, the deformation of the penetrating rings is particularly strong. The 
effective description, on the other hand, breaks down if the internal configurations of the rings in 
the monomer-resolved system differ significantly to the internal configurations in the system with 
only 2 ring polymers. Therefore, the anisotropic effective model should not be expected to be a 
quantitative description at the high densities in which this phase is formed. Agreement with the 
monomer-resolved model here is less satisfactory but the improvement over the isotropic model is 
still spectacular. 

For the rings with N — 100 we find that anisotropy does not play a key role any more, at least 
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Figure 9: The pair correlation function g(r) for a simulation of many ring polymers with N — 50 
monomers in the full monomer-resolved simulation (symbols), the anisotropic effective model 
(solid line) and the isotropic effective model (dashed line). The two plots show different reduced 
density p* ranges. 
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Figure 10: The pair correlation function g(r) for a simulation of many ring polymers with N = 100 
monomers in the full monomer-resolved simulation (symbols), the anisotropic effective model 
(solid line) and the isotropic effective model (dashed line). The two plots show different reduced 
density p* ranges. 
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not for the full model at the investigated densities. This had to be expected, as we could already see 
in figure [3] and [5] that anisotropy is less pronounced for larger ring sizes. As we saw in figure [6} c) 
both the isotropic and the anisotropic model give good results for g(r) up to p* ~ 2.5. In figure 
|T0} we see that for higher densities the inclusion of anisotropy does not yield results that are in 
better agreement with the full monomer-resolved simulations. The results in the isotropic effective 
model even seem to be in better agreement with the full model, which is attributed to multi-particle 
interactions that can change the configurations of the large and therefore more deformable rings 
significantly. The already small correlation between the directors, which is present in the dilute 
case, might therefore be even smaller at high densities. In the anisotropic effective model, we then 
overestimate the angular correlations between the directors and arrive at results that can be slightly 
worse than those of the isotropic model. Interestingly at p* — 20.0, which is the highest density 
investigated, the anisotropic model appears to crystallize. At this density we see the emergence 
of columns that are closed over the periodic boundary conditions and organize in a hexagonal 2D 
lattice structure. 

Finally, let us focus exclusively on orientational correlations. We define P(d f 11 ■ d f2 ^) as the 
probability density distribution for the scalar products between the directors d M 1 and d 2 of two 
ring polymers which are a distance r < 0.6D g o away from each other. In figure |TTj we present 
results for P(d^' • d 12 -’) for simulations in the monomer resolved and in the anisotropic effective 
model. If the directors were uncorrelated P(d^ • d f2) ) would be equal to 1. For low densities p* 
we obtain good agreement for all ring sizes investigated. Since we only look at the directional 
correlation of close by ring polymers, P(d 111 • d :2) ) can show strong anisotropic features even for 
p* —> 0. As expected the anisotropy in P(d 1 * • d i2 ^) is stronger for smaller rings. When the density 
is increased, the distribution always shifts towards parallel configurations in the effective model. 
This happens because less volume per ring is available for higher p* and by aligning parallel 
the rings occupy less space. Typically one observes the same trend in the monomer resolved 
simulation, only for the N = 100 rings we find more parallel rings for p* — 2.5 than for p* — 17.0. 
In contrast to the effective model the rings in the monomer-resolved simulation can deform and 
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Figure 11: P(d (1) -d* 2 )) is the probability density to find the scalar product d (1 ^ -d^ between the 
directors of two close by rings (r < 0.6D g o). Here we show P(d ( 1 ' • d (2 ^) for a simulation of many 
ring polymers in the full monomer-resolved simulation (symbols) and the anisotropic effective 
model (solid line), (a) N = 20; (b) N = 50; (c) N = 100. 


























their interaction with other rings can therefore be more isotropic at higher densities p*. This 
explains why for the large rings with N = 100 monomers, which deform more easily than the 
smaller rings, the correlation between the directors is much weaker than in the effective model and 
can even decrease with density. For N = 50 one can see that the number of orthogonal rings in the 
monomer-resolved model at high densities is significantly larger than in the effective simulation. 
As described in ref. 52 58 for N — 50 and p* > 12.8 one observes that oblate rings are interpenetrated 
by elongated prolate rings. Since the directors of the oblate and the interpenetrating prolate rings 
can be orthogonal to each other, one observes perpendicular directors for N = 50 even at the highest 
densities investigated. In the anisotropic effective model on the other hand, this interpenetration is 
disfavoured and we observe almost no orthogonal close-by rings at p* = 20 for N = 50. 

6 Truncation of the Expansion of the Anisotropic Potential 

Instead of working with a fully tabulated effective potential on a four-dimensional grid, it can 
be advantageous to use the analytical expansion on basis functions presented in the Appendix. 
Such expansions are truncated after some term, and here we shortly discuss the quality of such 
truncations for the problem at hand. To test the quality of the expansion of g (r,cos 0 \ .cos 62, <p) 
we also carried out Monte Carlo Simulations, where we used the effective potential associated 
to the expanded correlation function as the pair-interaction between our effective particles. We 
took the 14 coefficients C / 1 / 2 ,„(r) for which l\ .J2 < 4 into account and truncated the rest of the 
expansion. While g(r,cos0i,cos 02,9) can never be negative, the truncated expanded version of 
g can accidentally become smaller than zero. Wherever this happens g = 0 and V e ff = °o is used 
in the simulation. In figure [12] the pair correlation function g(r) obtained in this simulation is 
shown in comparison with the g(r) function, which we computed previously employing the full 
anisotropic effective potential. For both N — 20 and N — 50 we obtain reasonable results with 
the truncated effective interaction, given by only 14 expansion coefficients. For the full effective 
interaction, which we store on a 4D grid, we save 16 3 = 4096 entries for each value of r (see 
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Figure 12: The pair correlation function g(r) for a simulation of many ring polymers in the 
anisotropic effective model. For the dashed line we expanded the pair-correlation function g before 
computing the associated effective pair-potential. For the expansion we took the 14 coefficients 
for which l\, h < 4 into account. The solid line shows results of a simulation with the unexpanded 
effective pair interaction, (a) N — 20; (b) N = 50. 
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section 3.1). At intermediate densities the results obtained with the expanded effective interaction 
are a significant improvement with respect to the isotropic effective model. However, one has to be 
aware that for N — 20 the coefficients of higher order modes can still be quite high, especially for 
r between 0.2D g o and 0.7D a o. In figure [5]we see that the coefficient for the mode with {l\-h-m) = 
(0,2,0) can be larger than the coefficient of the isotropic expansion mode. The reason for the 
high contribution of higher order modes for N — 20 is of course the strongly peaked nature of 
g (r, cos 0i, cos 02, (p) for these small rings, which we can also observe in figure [3} The convergence 
of g (r, cos 0i, cos 02, (p) is poor for the N — 20 rings due to the strong anisotropy of their effective 
interaction. However our results show that the expansion modes up to l\,l 2 < 4 already capture 
the main features of the effective interaction. For N = 50 the degree of anisotropy is weaker and 
therefore the convergence of the expansion of g (r, cos 0i, cos 02 , (p) is better. 


7 Conclusions 

We have introduced a minimal anisotropic model to coarse-grain ring polymers with a finite bend¬ 
ing rigidity as soft, penetrable disks. For the shortest (N = 20) and the intermediate (N — 50) sized 
rings, this model represents a dramatic improvement over the isotropic coarse-graining, in which 
the relative orientations between the rings are all integrated upon and a radially symmetric interac¬ 
tion results instead. The approach is capable of distinguishing between the relative orientations at 
infinite dilution and it carries this distinction also to highly concentrated systems, where it repro¬ 
duces well the salient features of the structure as seen in the full, monomer-resolved simulations. 
Whereas this is valid more for N = 20 and N = 50, which have a contour length to persistence 
length ratio of N/s p ~ 2.7 and 6.7 respectively, some important features, such as the penetration of 
elongated rings in columns formed by oblate rings (found for N — 50), are suppressed or even lost 
in the effective description, as genuine many-body effects come into play. For the largest rings, 
N = 100, for which we obtain N/s p ~ 13.3, the contour length is much larger than the persistence 
length and they thus resemble more flexible objects. In this case the anisotropic potential at high 
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concentrations fails to describe the structural correlations. This indeed reflects the fact that such 


rings undergo, at high concentrations, conformational changes (shrinking, interpenetration) that 
are quite distinct from the assumptions that go into the anisotropic, soft disc-model, rendering it 
thereby very inaccurate. We therefore expect that our anisotropic model yields quantitative results 
over a broad density range, for systems of polymer rings with a contour to persistence length ratio 
ofiV/sp < 10. 

Our work provides, thus, an accurate and efficient general scheme for the coarse-graining of 
semiflexible ring polymers, as it allows for a very dramatic reduction of their degrees of freedom, 
while at the same time introducing a realistic class of systems for which anisotropic generalizations 
of the ultrasoft, penetrable effective interactions are physically meaningful. Future work will focus 
on the investigations of the structural and phase behavior of mixtures of stiff rings and of the 
dynamics of the structure formation in the same. 
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A Calculation of the Anisotropic Effective Potential 

As discussed in section[2j we wish to sample P(r,cos 0i,cos 02, 9 ) for a system of two ring poly¬ 
mers, in order to obtain a numerical expression for the anisotropic effective potential between 
them. We know that for large values of r, when the polymers can not interact with each other, 
P will correspond to the ideal case, Py. Therefore the interesting configurations for us are those 
values of r that result in overlaps between the ellipsoids of gyration. We use a biasing potential 
between the centers of mass of the two rings to restrict r to certain umbrella windows: 

v t,ii( r ) = ^ 2 • as) 
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With rj, kj we can tune respectively the location and width of the window for r in which con¬ 
figurations are sampled. We carry out simulations for a range of different rj and kj values and 
calculate histograms p[£ s (< 2) in the effective coordinates. Here, Q stands for (r.cos Q\ .cos 66 • <p) 
as a collective variable and thus denotes a bin in the effective coordinates, whereas j is the index 
of the respective biased simulation and thus determines kj and rj. The binning in the 4D space is 
identical for all biased simulations. We use the Self-Consistent Histogram Method by Ferrenberg 
and Swendsen 6465 to combine the different (<2), which results in an estimate P es tj(<2) for the 
histogram P (Q) of the unbiased system. The starting point of this method is that every simulation 
does in principle give an estimate for the histogram P(<2) of the unbiased simulation: 


= W“exp (PV<,{>JQ)) pW (Q). 


,U) 


(17) 


Here, denotes the bias potential in the j-th simulation, given by (16) and N {J > is a normaliza¬ 


tion factor, which can be expressed as 


ArM = £p( G )exp(-j3v“ (g)), (18) 

Q 

assuming that both P(<2) and Pj,-.| s (<2) are normalized. However, this estimate for P(<2) will only be 
useful for Q bins that have good statistics in the j-th simulation, which in our case means that the 
bins are at an r coordinate that is close to the rj value of the respective bias. Another problem with 
this expression is that in order to calculate N i j> we already need to know the sought-for quantity 
P(<2). To deal with the first problem, we combine the individual estimates obtained from each j-th 
simulation, to form an improved estimate: 


Pe S t(e) = £^ ) (G)p esU (Q). d9) 

j 

With c^'(Q) we can tune the weight of P e stj(<2) in the P (<2) estimate. We require Y.j c ^{Q) — 1- 
In bins where the j-th simulation has bad statistics we will choose c^> ( Q ) close to 0, such that the 
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Pest./ (2) estimate contributes only in bins where it is useful. The error of P est (Q) can be estimated 
via the Poisson distribution and it can be minimized via the following choice for the c2) (Q): 


exp(—/3 v[£ s (2))mWa^ 

E,exp(-/3(G))ilf(W*)' 


( 20 ) 


Here M ij) is the number of uncorrelated configurations sampled in the j-th simulation. With Eqs. 
we now arrive at an expression for an estimate of P(<2)- However, as N (j> depends on P(<2) 
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the expression (18) can only be evaluated if P(<2) is known in the first place. We can deal with 


this problem by using P es t(<2) for P(<2) i n the formula for N-^ (18) to obtain a self-consistency 
problem for P es t(<2), which can be solved iteratively. As a starting point for this iterative procedure 
an initial guess for P es t(<2) has to be provided. However, after many iterations the procedure 
is expected to converge to the same distribution, independent of the given initial condition. We 
started the iterative algorithm with an uniform distribution for P est (<2). 


B Expansion of g(r,cos 0i,cos 02,9) 

As the anisotropic effective potential V e ff (r, cos 0 \, cos On, (p) and the corresponding pair-correlation 
function g (r,cos 0\ .cos 0 2 , (p) depend on 4 variables, it is hard to visualize them. Nevertheless we 
can obtain a quantitative measure of the anisotropy in g and V e ff by carrying out an expansion of 
g: 


g (r, cos , cos 0 2 , (p) — £c„(r) f„ (cos 0i, cos 0 2 , <P) • (21) 

n 

Here, f, 7 are modes that depend on the angular degrees of freedom only and they form a complete 
basis for the angular dependence of g. 

To obtain a suitable set of basis functions f n for the dependence of g on the directors d (1,2 ) for 
a given vector r between the rings, we started with an expansion to a sum of products of spherical 
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harmonics: 


oo Z] OO l 2 

r (d< 1 ),di 2 ))= £ £ £ £ Q 1 ,m 1 ,/ 2 ,-» 2 y T 1 ( cos • fOy” 2 (cose 2 , n)■ oz) 

l\=Om\=—l\ h=0ni2=—l2 

Here 7 denotes g at a given r. By cos 0, and (p, the director d^ is represented in spherical coordi¬ 
nates. We use a reference frame where the connection vector r between the two rings points to the 
north-pole and therefore denotes the azimuthal angle around r. With Y! n we denote the spherical 
harmonics 66 . They fulfil the orthonormality relation f dQY"‘(cos 0, (p) Yjf (cos 0, (p) — 8i.i'8 m . m ' 
with d£l — d cos Od(p, where z denotes the complex conjugation of a complex number z. This 
allows us to compute the expansion coefficients via integration: 

c h,m u l 2 ,m 2 = J d£l\d(d. 2 Y™ 1 (cos 0 i, <jPi) Y ™ 2 (cos 0 2 , (pi) 7 (cos 0 i , <pi,cos 0 2 , (Pi)- (23) 

Due to the symmetries of 7 , we will now be able to compute or relate many of the coefficients 
and thus arrive at a reduced set of basis functions (f/, j 2: m} with which we can still represent 7 
exactly. We first use the continuous symmetry under rotations around the connection vector r: 


7 (cos 0 i,<pi,cos 0 2 ,<p 2 ) = 7 (cos 0 i, <pi — <p 2 , cos 02,0). 


(24) 


Using (23) and K/'^cos 0.(p) exp(wz<p) one can show that r :/,,/ 2 .„, 2 vanishes for in 1 ^ —m 2 


due to this symmetry. In the following, we enumerate the expansion modes with m = m\ — —m 2 . 
Next we use that d^ is equivalent to -d''-’ and therefore 7 fulfils the symmetry 7 ^d ! 11 .d 12 ' j = 
d (1) ,d (2 ^ = 7^d^\ — d (2 ^. Since T/"(d) = (—1 ) 7 U/” (—d), c / |mi / 2 „, 2 are zero if either l\ or 
h are odd. 

The monomer-resolved model is symmetric under a mirror transformation, which is therefore 
also a symmetry of the effective model. If we consider a state with <p 2 = 0 and mirror it by a plane 
spanned by r and d* 2 ^ we obtain a state with identical r, d (2 * and cos 0 i, while (pi changes sign. 
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Hence we obtain: 


/(cos 01, <Pi,cos 02, 92 ) = y(cos 01,9l - Cfh, COS 02, 0) = 


y(cos 0 i, 9 2 - 9 i,cos 0 2 ,O) = y(cos 0 i,<p 2 ,cos 0 2 , 9 i). 


(25) 


For the first and the last step we used (24). Therefore y is invariant under exchanging (p\ and 92 . 


We therefore have ci um j 2 _ m = c, u - m j 2 , m . 


Everything discussed so far also holds if we calculate the effective interaction between different 
rings, e.g. for rings with a different number of monomers. The final symmetry, which we will 
exploit now, only holds if we have identical rings. In this case we obtain an equivalent state if we 
swap the orientations of the two ring polymers: 


y(cos 0 i, <j?i, cos 0 2 , 92 ) = y(cos 02 ,<p 2 ,cos 0 i,<pi) = y(cos 0 2 , <pi,cos 01 , 92 ) 


(26) 


For the last transformation we used (25). Hence for identical rings y is also invariant under an 
exchange of 6\ and 02 and we therefore know that c / |m / 2 m — 

We now group basis functions if we a priori know that y has identical coefficients ci um j 2 m 
with respect to them. We sum the modes in each group and divide by a/#, where # is the num¬ 
ber of modes in the group. For identical rings # is at most 4, but can also be smaller if l\ — I 2 
or m = 0. For different rings # = 2 if in =4 () and 1 otherwise. In this way we obtain a new 
set of basis functions fi u i 2 , m • The indices of //, refer to the indices of one of the modes 
in the group from which fi u i 2j , n was constructed, with the additional constraint that m > 0 and 
h > h in the case of identical rings. As an example we can consider f2AA , which is set to 
(F 2 ' y-' +Y-' Yj + Yj y 2 ~ 1 + y 4 ~ 1 y 2 ') /V4 in the case of identical and (K 2 'T 4 1 +F 2 “ 1 T 4 1 ) /\/2 in 
the case of different rings. The new basis functions fi u i 2 , m °nly depend on (p = \(pi — 92 1 an d 
not on 91 , 92 separately. Thus 1 91 — 92 1 is precisely the 9 coordinate defined in ([ 6 ]), on which g 
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depends. Expanding the angular dependence of g with our new basis functions, we obtain 

min(/i,/ 2 ) 

g (r, COS 01, COS 02, p) = H, H H c h,l 2 ,m( r )fli,l2,m (cos 01, cos 0 2 , 9 ) .(27) 

h ={0,2,4,...} / 2 ={/ 1 ,/ 1 +2,/ 1 +4,...} m =0 

In the case of different rings, the sum over Z 2 does not start at l\ but at 0. Like also 

{f l\ ,/ 2 .m } fulfil the orthonormality relation 

Id£l 1 d£i.2 fi\,i 2 ,m (cos 0!,cos 0 2 , 9 ) fi^ m !(cos 01 ,cos 0 2 , 9 ) = Sij'Smrf. (28) 


Note that complex conjugation is not necessary as fi u i 2 . m are real functions in contrast to Y\ n Y ] 


h h 


Hence, we can obtain the coefficients c i u i^ m (r) via an integration analogous to (23): 



dQ.idQ .2 fi u i 2 , m (c os 0 i, cos 0 2 , 9 ) g(r,cos 0 ],cos 0 2 , 9 ). 


(29) 


To calculate c/ b / 2jm (r) we could do a numerical integration of g(r.cos 0\ .cos 0 2 . (p). A different 


approach for the calculation of c/ 1; / 2;W (r) is to express the integral in (29) as an average over con 


figurations of a system of 2 ring polymers for fixed r. We sample these configurations from the 
simulations with the bias potential Vbi as (/) given in (jT6]), which we carried out for calculating g on 
a 4D grid using umbrella sampling. As Vbi as ( r ) does not change the relative weight of configura¬ 
tions with identical r values, we can use the configurations that fall in a small window of r values 
to estimate averages over the angular degrees of freedom at some fixed r value. The average of the 
expansion modes //, j 2Jn over the configurations of the rings is related to C[ t ,i ljn (r) as follows: 


{fll,l 2 ,m)r ~ 


f dQ.idQ .2 g (r, cos 0 i, cos 0 2 , 9 ) f h (cos 0 i, cos 0 2 , 9 ) 


/ Jf2i<if2 2 g(r,cos 0i,cos 0 2 , 9) 


Cl u l 2 ,m{ r ) 

( 47 T) 2 g iso (r 


(30) 


Here we used (|9J) and (|29j) for the final step. g lso is easy to calculate numerically once we know 
g. The advantage of this approach with respect to a numerical evaluation of ([29]) is that we do not 
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need to introduce a grid in the angular degrees of freedom. In particular if we want to calculate 
c h,l 2 ,m( r ) coefficients for high l\ , I2 values the correctness of the result in the first approach depends 
sensitively on the number of grid points and also on the numerical method for carrying out the 
integration. Calculating averages on the other hand is straight-forward and we can easily estimate 
the statistical error using the standard deviation of block-averages. For these reasons we calculate 
the expansion coefficients c/ 1; / 2 , m (r) and estimate their errors following the latter approach. 

The larger the faster the fi u i 2jm functions can change when the director dW is varied. /o,o,o = 
(4/rj~ 1 is constant and therefore gives the isotropic contribution. By comparing co,o,o(t) with the 
size of the coefficients for (/ 1 , L.m) =4 (0,0,0) we can quantify the importance of anisotropy in the 
effective interaction. 

We find that with the 14 coefficients C/ 1; / vra (r) for which li,h < 4 we already get a reason¬ 
able approximation of g. This is true, even for the smallest rings investigated (TV = 20), where 
the anisotropy is most important. This fact allows us to store the essential information in the 
anisotropic potential with only a few functions of one variable, c/ t ; 2 m (r), instead of a 4D grid with 
a very large number of grid points. Using only the coefficient co.o.o(r) we recover g 1S0 (r) of the 
isotropic effective model. 
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